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Abstract. We present a study of the mixing properties of the simulated intra cluster medium, using tracers particles that are 
advected by the gas flow during the evolution of cosmic structures. Using a sample of seven galaxy clusters (with masses in 
the range of M ~ 2 - 3 • 10^^ Mq/H) simulated with a peak resolution of 25kpc/h up to the distance of two virial radii from 
their centers, we investigate the application of tracers to some important problems concerning the mixing of the ICM. The 
transport properties of the evolving ICM are studied through the analysis of pair dispersion statistics and mixing distributions. 
As an application, we focus on the transport of metals in the ICM. We adopt simple scenarios for the injection of metal tracers 
in the ICM, and find remarkable diff'erences of metallicity profiles in relaxed and merger systems also through the analysis of 
simulated emission from Doppler- shifted Fe XXIII lines. 



Key words, galaxies: clusters, general - methods: numerical - 
intergalactic medium - large-scale structure of Universe 



1. Introduction 

The mixing properties of the intra cluster medium (ICM) are 
still poorly understood and a number of open fields of research 
are tightly connected to this important topic: the 'heating cool- 
ing flows' problem (e.g. Bruggen & Kaiser 2002; Omma et 
al. 2004), the spreading of metals in the ICM (e.g. Schindler 
& Diaferio 2008; Borgani et al.2008), the stability of ther- 
mal profiles in the innermost region of galaxy clusters (e.g. 
Sharma et al.2009), the efi&ciency in the ram pressure strip- 
ping of accreting satellites and the emergence of cold fronts in 
galaxy clusters (e.g. Ascasibar & Markevitch 2006; Markevitch 
& Vikhlinin 2007) and many others. 

From the theoretical point of view, in the last few years 
meaningful progress has been made in understanding the con- 
vective properties of the ICM. Focusing on the role played 
by instabilities in a magnetized, low-collisional plasma (e.g. 
Parrish & Stone 2005; Quataert 2008; Ruszkowski & Oh 2009), 
a new picture of the ICM stability was presented, which drasti- 
cally alters the "classic" picture provided by the Schwarzschild 
criterion, in which stability is ensured by dlog(S)/dr > 
(Schwarzschild 1959), where S is the specific entropy. 

These approaches cannot take into account all sources of 
mixing motions in galaxy clusters though, which are due to 
large scale accretions of matter and turbulent motions. Indeed 
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there is clear evidence nowadays that a sizable amount of tur- 
bulent motions may be present in the ICM. Observations sug- 
gest large scale subsonic turbulent motions in the range of 
~ lOOkpc - IMpc (e.g. Schuecker et al.2004; Henry et al.2004; 
Churazov et al.2004). In addition, also studies of Faraday 
Rotation allow a complementary approach and suggest that the 
ICM magnetic field is tangled over a broad range of scales (e.g. 
Murgia et al.2004; Vogt & Ensslin 2005; Kuchar & Ensslin 
2009; Bonafede et al.2010); also, the diflTuse radio emission 
detected in a fraction of merging clusters may provide indirect 
evidence for turbulence in the ICM (e.g. Brunetti et al. 2008). 

Obtaining a self-consistent picture of the evolving ICM, 
where large scale and small scale mixing properties of galaxy 
clusters are followed during the whole cluster evolution is chal- 
lenging, and in this respect cosmological high resolution nu- 
merical simulations may provide additional and valuable infor- 
mation. 

At present, two main numerical methods are massively 
applied to cosmological numerical simulations: Lagrangian 
methods, which sample both the Dark Matter and the gas 
properties using unstructured point-like fluid elements, usu- 
ally regarded as particles (e.g. smoothed particles hydrody- 
namic codes such as GADGET, Springel 2005, or HYDRA, 
Couchman et al.l995) and Eulerian methods, which recon- 
struct the gas properties with a discrete sampling of the space 
using meshes of various resolution (e.g. cosmological TVD 
code, Ryu et all 993; ENZO, Norman et al.2007; FLASH, 
Fryxell et al.2000; RAMSES, Teyssier 2002 etc) and model the 
Dark Matter properties with a particle mesh approach. 

The typical advantages and drawbacks of the two ap- 
proaches have been extensively investigated. Lagrangian meth- 
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ods allow one to obtain a high spatial resolution where matter 
concentrates, but in low density regions the sampling is rather 
poor. By construction, this class of methods allows one to fol- 
low in a natural way the advection of a single parcel of (gas 
or DM) matter during the whole simulation and its dynamical 
evolution. 

Eulerian methods have a spatial resolution fixed to the size 
of the cell in the computational mesh, and this often is inade- 
quate to follow the details of cosmological and physical pro- 
cesses of interest in the innermost region of collapsed objects. 
The adaptive mesh refinement technique (AMR) overcomes 
this problem by increasing the spatial resolution of the Eulerian 
mesh in regions of particular interest (e.g. Berger & Oliger 
1984; Berger & Colella 1989; Kravtsov et al.l997). A further 
important feature of Eulerian methods is that they are strictly 
conservative for the total fluid energy, and therefore they ac- 
curately follow the Rankine-Hugoniot relations for shocks (a 
property shared by shock capturing methods, also known as 
Godunov schemes, Godunov 1959). This is particular impor- 
tant in cosmological simulations, where cosmic shocks play an 
important role in setting both thermal and non thermal prop- 
erties of the ICM (e.g. Ryu et al.2003; Pfrommer et al.2006; 
Vazza et al.2009). 

In the last few years a number of works highlighted and 
investigated the causes which lead to the main difl'erences be- 
tween these numerical approaches (Agertz et al.2006; Tasker 
et al.2008; Wadsley et al. 2008; Mitchefl et al.2009; Robertson 
et al 2009; Springel 2009). The adoption of artificial viscosity 
by SPH and the limited spatial resolution for AMR, were the 
reasons for the more significant diff'erences between the two ap- 
proaches; if however both sources of error are tuned appropri- 
ately (e.g. by reducing the action of artificial viscosity outside 
shocks in SPH, or by increasing the mesh resolution in AMR 
codes), the two approaches produce consistent results. 

High resolution AMR simulations (such as the ENZO sim- 
ulations presented in this paper) can provide an accurate repre- 
sentation of the cosmic gas dynamics in galaxy clusters, achiev- 
ing a very broad dynamical range. Recent works showed that 
opportune refinements (e.g. lapichino & Niemeyer 2008, Vazza 
et al. 2009; Maier et al.2009; Zhu, Feng & Fang 2010) allow for 
efl&ciently studying the onset and evolution of chaotic motions 
in the ICM. 

Yet, some interesting information of Lagrangian nature 
cannot be followed by these methods. One example is the 
distribution of metals in the ICM due to diff'usion and turbu- 
lent mixing which cannot be followed, unless the conservation 
equations are self-consistently implemented for every metal 
species, with a sizable increase in algorithmic complexity and 
computational expense. However, it is possible to treat metals 
in a post-processing stage by following their spatial evolution 
through that of mass-less particles (tracers), which move along 
with the baryon gas. 

An other interesting case is the advection of cosmic rays 
(CR) particles in the ICM. Shocks, turbulence, collisions be- 
tween high energy hadrons and AGN/supernovae activities are 
expected to inject relativistic particles in the ICM over cosmo- 
logical time-scales (e.g. Brunetti & Lazarian 2007), whose in- 



terplay with diff'use ~ \iG ICM magnetic fields is responsible 
for the diff'use non-thermal radio emissions observed in galaxy 
clusters (e.g. Ferrari et al.2008 for a recent review). If the en- 
ergy stored in CR is much smaller than the thermal ICM (as 
shown by recent observations of the central regions of clusters, 
e.g. Reimer 2004; Aharonian et al.2008; Brunetti et al.2007; 
THE MAGIC collaboration et al.2009), the spatial evolution of 
CR can be studied through that of passive tracers. This applies 
as long as the eff'ect of CR diff'usion is negligible and the evo- 
lution of CR particles can be regarded as the simple advection 
problem of tracers in the evolving ICM. 

The objective of the present work is to investigate the mix- 
ing properties of the ICM with an appropriate numerical recipe. 
For this purpose high algorithmic accuracy and high spatial 
and mass resolution are needed. These requirements are well 
provided by ENZO, which is a high-order and cosmological 
AMR Eulerian code (e.g. Norman et al. 2007). Here we present 
ENZO simulations with two important customizations: first, an 
additional refinement criterion is added to increase the spatial 
resolution on propagating shocks (Vazza et al.2009); second, 
since a pure Eulerian method cannot provide the details on the 
behavior of each specific fluid element (e.g. its trajectory and 
velocity), we explicitly follow the the advection of a large num- 
ber of Lagrangian passive (mass-less) tracer particles, which 
move according to the 3-D velocity field of the gas. 

The paper is organized as follows: in Sectl2]we present the 
details of the simulation runs for this project, and the method 
to inject and follow tracers. In Sectl3]we discuss the result of 
several convergence tests to assess the reliability of tracers sim- 
ulation with diff'erent possible setups. In Sect|4]we present as- 
trophysical results of tracers, studying in particular the average 
transport properties of tracers in all simulated clusters by dis- 
cussing the evolution of the pair dispersion statistics (Sect j4.1l) : 
the morphology and the evolution of radial mixing in merg- 
ing and non merging clusters through the evolution of tracers 
initially located in spherical shells (Sect j4.2l) : the spreading of 
metal tracers in the ICM, studying toy models of metal injec- 
tions from galaxies (Sect j4.3l) : the simulated emissivity from 
the broadened Fe XXIII line from metal tracers and its depen- 
dence on the the dynamical history of clusters (Sect l4.4l) . Our 
conclusions are given in SectJS] 



2. Numerical methods 

2A. The ENZO code 

ENZO is an AMR cosmological hybrid code originally written 
by G. Bryan and M. Norman (Bryan & Norman 1997; Norman 
et al.2007). 

It uses a particle-mesh N-body method (PM) to follow the 
dynamics of the collision-less Dark Matter (DM) component 
(Hockney & Eastwood 1981). The DM component is coupled 
to the baryonic matter (gas), via gravitational forces, calculated 
from the total mass distribution (DM-hgas) solving the Poisson 
equation with a EFT based approach. 
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Fig. 1. An example of the advection of tracers within a forming galaxy cluster. Contours: projected gas density across a volume 
of side 10 Mpc/h. Blue points: projected positions ofN^ 10^ at z = 21, z = 3 and z = 0. 



Table 1. Main characteristics of the simulated clusters. Column 
1: identification number. Columns 2: total mass (Mdm + ^gas) 
inside the virial radius. Columns 3: virial radius. Column 4: dy- 
namical state at z = (qualitative classification); MM=major 
merger for z < 1; mm=only minor mergers for z < 1, 
rr= visually relaxed at z ~ 0. Cluster Hid is the same as HI, but 
re- simulated with mesh refinement on gas/DM density only. 



ID 


Mw.[1014Mo//z] 


R^kpc/h] 


dynamical state 


HI 


3.10 


1890 


mm 


Hid 


3.10 


1890 


mm 


H2 


3.05 


1810 


MM 


H3 


2.95 


1710 


rr 


H4 


2.63 


1738 


rr 


H5 


2.41 


1703 


MM 


H6 


2.32 


1593 


MM 


H7 


2.14 


1410 


mm 



The gas component is described as a perfect fluid, and 
its dynamics are calculated by solving the conservation equa- 
tions of mass, energy and momentum over a computational 
mesh with an Eulerian solver based on the Piecewise Parabolic 
Method (PPM, Woodward & Colella, 1984). The PPM algo- 
rithm belongs to a class of schemata in which shock waves 
are accurately described by building into the numerical method 
the calculation of the propagation and interaction of non-linear 
waves with a directionally split approach. It is a higher or- 
der extension of Godunov's shock capturing method (Godunov 
1959). It is at least second-order accurate in space (up to 
the fourth-order in 1-D, for smooth flows and small time- 
steps) and second-order accurate in time. The PPM method 
requires no artificial viscosity, which leads to an optimal treat- 
ment of energy conversion processes, to the minimization of 
errors due to the finite size of the cells of the grid and to a spa- 
tial resolution close to the nominal one. The basic PPM tech- 
nique has been modified to include the gravitational interac- 
tion and the expansion of the Universe (e.g. Bryan et al.l995). 
In ENZO cosmological simulations, several comparison tests 
(e.g. O'Shea et al. 2005; Tasker et al. 2008)) provided evidence 
of the better performance of the PPM method implemented. 



compared to the alternative method of ZEUS artificial viscos- 
ity (which is also implemented in ENZO). 

2.2. Cluster simulations 

For the simulations presented here, we assumed a "concor- 
dance" ACDM cosmology with the parameters Qq = 1-0, 
QfiM = 0.0441, Dz)M = 0.2139, Qa = 0.742, the Hubble pa- 
rameter h = 0.72 and a normalization for the primordial density 
power spectrum erg =0.8. 

The clusters considered in this paper were extracted from 
independent cosmic volumes, each with the size of 109Mpc/h, 
obtained with root grids of 256^ cells and with DM mass reso- 
lution of ntdm ~ 5.4- lO^M© /h. We additionally sub-sampled cu- 
bic volumes of side 54.5Mpc/h with a second 256^ grid around 
the region of the formation of our clusters, achieving an eff'ec- 
tive root grid resolution of ~ 2l0kpc and an eff'ective DM mass 
resolution of mdm ~ 6.7 • 10^ Mq/H. For every cluster run, we 
identified cubic regions of the size of ~ 47?v/r (where Ryir is the 
virial radius of clusters at z=0), and allowed for three additional 
levels of mesh refinement, achieving a peak spatial resolution 
ofA^25kpc/h^ 

The mesh refinement was triggered by gas or DM over- 
density criteria from the beginning of all simulations (z = 30). 
From z = 2 an additional refinement criterion based on 1-D ve- 
locity jumps was switched on (Vazza et al.2009). This second 
AMR criterion is designed to capture shocks and intense turbu- 
lent motions in the ICM out to the external cluster regions. The 
grid was refined of a factor two whenever one of the following 
conditions was matched: 

- either the gas or the DM density at a given cell exceeded 
^p/Pmean > ^p, whcrc p is the gas (DM) density within the 
cell at the 1-AMR level, and Pmean is the mean gas (DM) 
density at the root grid level (in the AMR region); 

- the 1-D velocity jump across a patch of three cells (normal- 
ized to the minimum velocity modulus in the same patch) 
was larger than a threshold \Svi\/min(vi) > Sy, where i = x,y 
or z. 

^ Below we will refer to this sub-volume as to the 'AMR region'. 
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Fig. 3. Projected paths of four random tracers, evolved from z = 30toz = with the NGP method (red lines) and the CIC method 
(blue lines); the isocontours show the projected gas density at z = 0. The side of the image is lOMpc. 
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Fig. 2. Radial position of tracers in a simulated 2-D rotating 
disk for different time interpolation schemes: NGP with up- 
dates every time step of the simulation (solid lines), every two 
time steps (dotted), every 4 time steps (dashed) or with an im- 
plicit second order scheme (solid black). From top to bottom, 
results are shown for increasing mesh resolution; the horizontal 
red line shows the ideal results. 

The adopted threshold values to trigger AMR were Sp = 3 
and Sy = 3. The second AMR criterion was switched on only 

^ We note that this threshold values are slightly smaller than in other 
AMR simulations, e.g. Nagai et al.2007; lapichino & Niemeyer 2008. 



from z = 2, since this produced only a negligible difference 
in the cluster simulations compared to its application starting 
from larger redshifts (Vazza et al.2009). 

At z = 0, the number of cells refined up to the peak resolu- 
tion in our runs corresponded to a ~ 20 - 40 per cent of the total 
volume of the AMR region (A/^ ~ 10^ - 10^ cells). Compared 
to standard cluster simulations, where mesh refinement is trig- 
gered by gas/DM over-density, the dynamical range available 
for chaotic motions within R^ir is much larger, since turbulent 
motions with significant 1-D jumps in velocity are not artifi- 
cially damped by the under- sampling effects due to poor res- 
olution, which would otherwise arise if they are not following 
large enough over-densities (of gas or DM) in the cluster vol- 
ume. 

Table [T] summarizes the main properties of the galaxy clus- 
ters simulated in this work. The last column in the Table re- 
ports a broad classification of their dynamical state, based on 
the presence of merger events in the range < z < 1 . All runs 
presented here are non-radiative and no treatment of the reheat- 
ing background due to AGN and/or stars is considered. 

In order to understand the effect played by our refinement 
strategy on the properties of simulated tracers, we addition- 
ally re-simulated the most massive cluster of our sample using 
the standard refinement criterion (triggered by gas/DM over- 
density) starting from z = 30 (run Hid). A more general com- 
parison on the effect played by the refinement criteria on the 
Eulerian properties of the simulated clusters can be found in 
Vazza et al. (2009). 

2.3. Tracer simulations 

The transport and diffusion properties of turbulent media are 
naturally addressed in the Lagrangian viewpoint, and many 
highly sophisticated numerical techniques haven been devel- 
oped and tested over years (e.g. Toschi & Bodenschatz 2009 
and references therein). Usually, the statistical properties of 
fluid particles transported by a 3-D fully developed turbulent 
flow are investigated by means of direct numerical hydro sim- 
ulations and by following the advection of mass-less (or in- 
ertial) particles, whose 3-D position is timely updated by as- 
suming the motion of particle tracked by the flow (e.g. Maxey 
& Riley, 1983). In the last few years, a number of interesting 
works have investigated in detail the behavior of tracers in as- 
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Fig. 6. Black points: projected position for three different generation of tracers (generation '0' is generated at z = 30, gen. '3' at 
z ^ 6, and generation '9' at z ^ 0.1) for a cluster run at z = 0. Color maps: slice through the center of the cluster, showing gas 
density (top) and gas temperature (bottom). The side of the images is lOMpc. 



trophysical turbulent flows by taking into account the effects 
played by inertia effects (e.g. Bee et al. 2006; Calzavarini et 
al. 2008), intermittency (e.g. Biferale et al. 2004), gravity (e.g. 
Bosse et al. 2006) and magnetic fields (e.g. Arzner et al. 2006). 

In the framework of galaxy clusters study with cosmolog- 
ical numerical simulations, to our knowledge no study of the 
Lagrangian properties using tracers particles have been per- 
formed so far. In this work, we present a method to inject and 
follow tracer particles in cosmological adaptive mesh refine- 
ment simulations with the ENZO code. A qualitative similar 
analysis was recently presented by Federrath et al. (2008), who 



used ENZO PPM simulations to study the turbulent transport 
in the Interstellar Medium. 

We followed the kinematics of passive tracers by inject- 
ing and evolving Lagrangian mass-less points in the AMR re- 
gion of each cluster and by updating positions according to the 
underlying gas velocity field. Since the volume considered for 
tracer studies is always centered on a forming massive galaxy 
cluster, the bulk of tracers is progressively advected towards the 
center of the computational domain; in the (rare) case of tracers 
exiting from the considered volume, we simply removed them 
from the computation. 
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Fig. 4. Profiles of number distribution of tracers in run H5 at 
z=0, for tracers injected at z = 30 (blue) and tracers generated 
at z = 2 (red) . The continuous lines are for the CIC interpo- 
lation, the dotted lines are for the NGP interpolation and the 
dashed lines are for the CIC interpolation with time integration 
with coarse time steps. 



Fig. 5. Profiles of number distribution of tracers in run H7 
at z=0 (blue lines) and at z=1.0 (black lines). Tracers have 
been evolved using with NGP method with updates every two 
time steps (in solid) and with an implicit second order scheme 
(dashed). 



In order to better sample the entire cluster volume during 
time, we injected diff'erent populations of tracers in the AMR 
region. For the remainder of the paper, we will refer to tracer 
"generation" as to the injection of mass-less tracer particles in 
the AMR region at a given cosmic epoch. When new tracers 
were generated, they were placed at the center of cells and were 
uniformly distributed in the grid. If diff'erent mesh resolutions 
were present, as in the case we analyze here, the initial sam- 
pling was always taken from the most refined mesh. 

In our fiducial setup, tracers were moved by updating their 
comoving coordinates through: 

x(t-\- At) = x(t)-\-uAt. (1) 

The time step used to evolve the tracer positions is the time 
step of the simulation. In non-radiative simulations the time 
step in ENZO is computed as the minimum between three time 
steps related to three typical velocities: the maximum velocity 
module among Dark Matter particles, vdm, the maximum gas 
velocity in the grid, Vgrid, and the the expansion velocity of the 
Universe at each time, vh- The actual time step in ENZO is 
thus: 

= nc/max\vDM, Vgrid. v//|; (2) 

where nc (the Courant safety number) is customary set to 
0.4. The above choice ensures that no signal can travel for more 
than a half-cell during a time step (at all AMR level) and it also 
ensures that updating the tracer velocities every two time steps 
is accurate enough. 



In principle, more accurate time integration schemes can be 
adopted to better preserve stability, like the Runge-Kutta or the 
leapfrog scheme (e.g. Hockney & Eastwood, 1981); in Sect l3.ll 
we specifically present a comparison between diff'erent strate- 
gies for the time-integration of tracer positions. 

The tracer behavior was calculated through a post- 
processing procedure using the outputs of ENZO simulations 
as input for evolving the tracers positions in time. In principle 
the whole procedure can be incorporated in ENZO as an addi- 
tional run-time routine, but we found it much more convenient 
to follow a post-processing strategy, even if this led to a general 
overhead in terms of data storage and management. Indeed, the 
post-processing approach first allowed us to perform a large 
number of numerical tests to find the best setup for the gen- 
eration of tracers and for the computation of their evolution. 
Second, most of the physical mechanisms that we planned to 
simulate with tracers were expected to have negligible dynam- 
ical feedback on the surrounding baryonic matter (e.g. metal 
enrichment and spreading in the ICM, cosmic rays) and thus 
the same high-resolution cluster simulation can be used for a 
large number of studies. Finally, some of the physical mecha- 
nisms that can be studied with tracers, like the injection of CR 
at shocks, are still poorly understood theoretically have a still 
incomplete theoretical understanding and the iterative applica- 
tion of tracers may help to explore a wide number of diff'erent 
assumptions. 

As an example of a simple tracer run, we show in Fig[T]the 
projected positions for a single injection (generated at z = 30) 
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Fig. 7. Number density profile at z=0 for ten diff'erent gener- 
ations of tracers (one every ^ IGyr), the black curve is for 
tracers generated at z = 30, while the orange one is for tracers 
generated at z ~ 0.1. 

of N ^ 10^ tracers at three diff'erent epochs in a cubic sub- 
volume of the size of lOMpc/h. 

Section [3] presents numerical tests to decide the optimal 
setup of our tracer algorithm. In Sect j3.1l we compare the use 
of diff'erent interpolation scheme to assign tracers' velocities; 
in Sect l3.2l we discuss the eff'ect of diff'erent injection epochs; 
in Sect j3.3l we investigate the role played by the spatial sam- 
pling strategy; in Sect l3.4l we report on the suitable number of 
tracers that should be adopted to avoid undersampling of the 
Eulerian grid. 

3. Numerical tests 

3.1. What velocity for tracers? 

The numerical convergence of our method was tested by inves- 
tigating two approaches for the spatial interpolation to assign 
velocity to each tracer particle. 

The Nearest Grid Cell (NGP) scheme allots a tracer the ve- 
locity of the gas on the closest grid point. With the the Cloud 
in Cell (CIC) interpolation scheme though, the 3-D velocity at 
the tracer position is calculated as 

ijk 

where Vijk is the gas velocity and dijk is the linear distance 
between the center of the cells and the tracer, in units of the cell 
size (the sum is performed over the nearest eight cells). 

We tested the time interpolation procedure by implement- 
ing an implicit second-order method with a NGP assignment 
for the tracer positions: 

x(t + At) = x(t) -h (m -h u')At/2, (4) 



where m' is the 3D velocity field at the time t -h At. 

First, we performed an idealized 2-D test which compares 
the above interpolation schemes for a rigid rotating disk (FigIS 
for increasing resolutions of the underlying grid distribution 
(64^, 128^ and 256^). Tracers are initially located at a fixed ra- 
dius from the center of the disk and their motion is followed for 
1500 time steps (which is similar to the total number of time 
steps in typical clusters runs) and the time step At is timely 
computed with the same prescription from the Courant condi- 
tion of our ENZO simulations. For the NGP and CIC schemes 
we also tested the possibility of updating the tracer velocity ev- 
ery two or every four time steps, because this would reduce the 
computational cost in cluster simulations. For a perfect inter- 
polation procedure, tracers are expected to move in a perfectly 
circular way at the initial radius of injection. Figure [2] shows 
the time evolution of the radial position for a tracer in the ro- 
tating disk, placed in all cases at the initial distance of R/4 (R 
is the disk radius) from the rotation center. In this simple setup, 
the CIC and NGP schemes produce identical results and the 
diff'erent lines are superimposed for all choices of A^. At all 
resolutions, the adoption of a second order implicit method for 
the time interpolation leads to convergent results, with a rapid 
oscillations of < 2 cells around the radius of injection. The 
CIC and the NGP methods led also to similar radial oscilla- 
tions on larger periods if the tracer position was updated every 
At or 2At time steps, while they show an increasing radial dis- 
persion for the choice of 4 At at all grid resolutions. The above 
findings have been tested and confirmed even in the case of 
an angular velocity that varies with time (which is reasonable 
as the Courant condition is designed to preserve time accuracy 
against any change on the dynamical time of the system), and 
for diff'erent injection radii. Since the typical grid size corre- 
sponding to the AMR volume of our clusters is N ^ 200 - 300 
cells (at the highest resolution), this test suggests that updating 
tracer positions every < 2 time steps either with a CIC or with 
a NGP interpolation may be a viable choice for cluster simula- 
tion, where large rotation patterns may be present. 

To investigate this issue further, we simulated with the CIC 
and the NGP scheme the advection of tracers in a cluster run 
by using Eq[TJ Figure [3] shows the evolution of the projected 
positions of four tracers in a cluster simulation from z = 30 to 
z = 0. The color coding shows the path according to the two 
integration procedures: diff'erences in the positions are usually 
< 5 cells after ~ 150 - 200 integration time steps. 

Figure |4] shows the radial number density profile of trac- 
ers (i.e. mean number of tracers for unit of volume, for shells 
around the cluster center) at z = 0, adopting N ^ 10^ trac- 
ers that have been moved using both methods and Eq[T] for 
two diff'erent generation epochs, z = 30 (blue lines) and for 
z = 2 (red lines). We also show the distributions for an addi- 
tional run, which makes use of the NGP with positions updated 
every four time steps of the simulation. Statistically, the three 
methods provide consistent results, with no systematic diff'er- 
ences. Finally, in Fig {5] we show the results or our test of our 
fiducial NGP interpolation with updates every 2A^ against the 
implicit second-order interpolation scheme of Eq|4]for the ad- 
vection of N ^ 10^ tracers injected at z = 30 for the smallest 
cluster of the sample. The agreement at z = is better than 




Fig. 8. Tracer number profiles at z = for three diff'erent epochs of injection (zinj =10, Zinj = 2 and Zinj = 1); the dotted black 
lines show the distribution of tracers initially located in gas clumps, the dashed blue lines show the distributions of tracers initially 
located in a spatially uniform way, while the solid green lines show the average of the two; the solid red lines show the real gas 
density profile of the cluster at z = (as measured from the Eulerian cells). All runs have the same number of tracers; the gas 
density profile was re-normalized to the total number of tracers inside R^ir. 



~ 5 percent at all radii, while for z = 1.0 there is a larger 
scatter of up to a factor ~ 2 within lOOkpc from the cluster 
center. Given the considerably larger amount of data that the 
implicit second order method requires (since one must simulta- 
neously consider the 3-D velocity field at two close time steps), 
and given the relatively small statistical diff'erence compared to 
the other more simple schemes from now on we will adopt the 
NGP method, which seems to provide a sharper reconstruction 
of the gas velocity field close to shocks, and we will integrate 
the tracer positions every two time steps of the simulation. As 
a cross check, we report in Sect j4.1l an additional test to show 
that the results obtained with second order integration and sim- 
ple Euler step are very similar 0. 

3.2. The injection time. 

Matter accreted by clusters may present diff'erent dynami- 
cal properties according to the cosmic epoch of accretion. 
Consequently diff'erent generations of tracers will be spatially 
distributed depending on their injection epoch and on cluster 
dynamical evolution. 

For example, Figl6] shows the projected positions at z = 
of three diff'erent generations of tracers in a cluster run. Tracers 
were injected at z = 30, z = 6 and z = 0.1), with an initial 
spatial sampling of one tracer every four cells (in 1-D) of the 
grid at the maximum refinement level. 

Generations injected at earlier epochs are more concen- 
trated because of the cluster formation process. This is also 
evident from Fig|71 where we show the radial number density 
profile of tracers at z = for the ten different generations in the 
same cluster run. Diff'erences at the level of ~ 3 - 5 in the num- 
ber density profile are found within < lOOkpc from the cluster 
center, although all profiles have similar shapes, and a trend to 
present flatter distributions is found for the latest generations. 

^ Similar results on the convergence between Euler step and higher 
order time interpolation for ENZO simulations of the inter stellar 
medium were reported also in Federrath et al. (2008). 
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Fig. 9. Profiles of number distributions (bottom curves) and for 
average density (in arbitrary units, top curves) for tracers in a 
cluster run as a function of the total number of tracers. The 
black dashed line shows the volume weighted gas density pro- 
file measured in the Eulerian grid. 

3.3. The sampling strategy. 

Tracers basically represent a way to sample the Lagrangian 
evolution of the physical quantities associated to a fluid de- 
scribed with an Eulerian method, and the choice of the sam- 
pling strategy to place tracers in the mesh cells depends on the 
physical quantity of interest. For example, if the evolution of 
gas matter accretion is the quantity of interest, tracers should 
sample the Eulerian space in a density-weighted way rather 
than in a spatial uniform manner. In fact, a regular volume sam- 
pling of the grid would place most of tracers in smooth and 
under-dense cosmic regions, while a density-weighted sam- 
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Fig. 10. Phase distribution for four tracers runs, with decreasing sampling of the underlying Eulerian grid: from left to right, 
= 10^, = 10^ = 5 • 10^ and - 10^ tracers; the number profiles were normalized to the same total number for the sake 
of comparison. The isocontours show the phase diagram for the underlying gas distribution, with sqrt-spacing. 



pling would place most of them inside the self-bound gas 
clumps. For this reason the adoption of the initial sampling 
strategy (volume- weighted, gas density-weighted, DM density- 
weighted, etc) must be defined according to the investigated 
physical process. 

Figure IH] shows the results of a test where tracers were in- 
jected at diff'erent epochs (z = 10, z = 2 and z = I) following 
two diff'erent approaches: a) tracers were placed only inside the 
virial volume of the 50 most massive halos in the AMR region 
at the three diff'erent redshifts, with a number profile equal to 
that of the gas in each clump 0; b) tracers were placed with 
a regular spatial sampling of the grid, as in the previous sec- 
tions. Although this represents just a crude test to study the en- 
vironmental dependence of tracers, the trends found are clear. 
Tracers injected in clumps produce more concentrated number 
density profiles at z = 0, while tracers injected in the smooth 
gas component have lower densities in the core at z = and 
flatter profiles outside ~ O.SRyir. This tendency is increasingly 
clear as the redshift of injection decreases. For tracers injected 
at z = 10, the agreement of the average profile with the gas 
mass density distribution is generally within < 50 per cent. 
With a more accurate sampling of the gas density field at the 
epoch of generation, the gas density profile at z = can be 
closely reproduced by tracers. 

3.4. How many tracers? 

The choice of a suitable number of tracers is important to find 
the best compromise between computational time and accurate 
sampling of the underlying gas distribution. 

We performed convergence tests on the number of tracers 
with ten generations, with a number of tracers N = 10^, N = 
10^, A/^ = 5 • 10"^ and N = lO'^ for each generation. The time- step 
is kept the same in all tests. 

In Fig|9l we plot the radial number density profile of tracers 
and the average gas density profile at tracers positions. Tracers 
provide a good statistical sampling of the real gas density pro- 
file, with a scatter of less than < 10 per cent inside the virial 

^ Halo positions and virial parameters were computed according to 
a halo finder algorithm working with spherical gas+DM matter over- 
density, specifically designed for Eulerian simulations (e.g. Gheller, 
Pantano & Moscardini 1998; Vazza, Brunetti & Gheller 2009) 



radius. Also the number profile of tracers (when normalized to 
the same total value) shows convergent results, with the typical 
scatter of less than a factor ~ 2 within 500kpc from the cluster 
centers when the cases N = 10^ and A^ = 10^ are compared. 

Figure[TOl shows the density vs temperature phase diagrams 
for the four samplings discussed above(as colored points), 
compared with the phase diagram from the underlying distribu- 
tion of gas within the AMR region. The comparison shows that 
even if the tracer distributions agree well within the main halos, 
larger diff'erences due to poor sampling in low density regions 
emerge in runs with a sampling worse than one tracer every 
-100 gas cells. We therefore consider that, in order to have 
a fairly good representation of the cluster regions and also of 
the accretion pattern outside clusters, an initial sampling with 
at least 1 tracer every ~ 70 - 80kpc/h in 1-D (i.e. A^ ~ 10^ 
tracers in the AMR region) for every generation is necessary. 



4. Results 

4.1. Tracers dispersion statistics 

The presence of complex and turbulent velocity fluctuations 
in the simulated ICM has been established in a number of 
works which were performed with very difl'erent numerical 
techniques (e.g. Norman & Bryan 1999; Dolag et al.2005; 
Vazza et al.2006; lapichino & Niemeyer 2008; Lau et al.2008; 
Vazza et al.2009; Xu et al.2009). Using the same AMR tech- 
nique applied in this work, Vazza et al.(2009) were able to mea- 
sure the 3-D velocity power spectrum of a simulated galaxy 
cluster over a range of ~ 100 in spatial scales. Although cos- 
mological simulations provide a simplified view of the ICM, 
the presence of subsonic turbulent velocity fields in clusters is 
in line with existing observations (e.g. Schuecker et al.2004; 
Henry et al.2004; Sanders et al. 2009). 

Complex motions in the ICM aff'ect the mixing and trans- 
port processes in the baryon gas. In Fig[TT]we show the evolu- 
tion of the projected position of three random couples of trac- 
ers, injected at z = 30 with an initial separation of 36kpc. The 
tracers follow a rather laminar path in the first stage of their 
evolution, while their paths become tangled and fairly diff'erent 
when they enter the collapsing region of the forming cluster. 
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Fig. 11. Projected paths of three couples of tracers evolved from z = 30 o z = from an initial separation of 1 cell=36kpc, for 
run H7. Isocontours show the projected gas density at z = 0. The side of the image is 3.8Mpc 
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Fig. 13. Left .'Time evolution of the pair of tracers dispersion for A/^ = 2 • 10^ tracers in all runs for initial separations of 72kpc (two 
cells). The additional thick dashed line shows the slope for the Richardson law, R(t) ~ P^^. Right: 3-D velocity power spectrum 
of the galaxy clusters in the sample. The spatial scale k was normalized to the scale corresponding to the virial radius of every 
cluster ^0 ^ InjRyir. The additional dashed lines show the -5/3 scaling to guide the eye. Colors are the same as in the left panel. 



As a guide reference to understand the motions of tracers 
within the complex gas/DM velocity field in galaxy clusters 
one may use studies customarily performed for isotropic and 
incompressible turbulence (e.g. Bee et al.2006; Bee et al.2009 
and references therein). In this case the tracer pair dispersion 
statistics (i.e. P{t) =< \Xi{t) - X2{t)\ >, where Xi{t) and Xiit) 
are the position of couples of tracers at any time) are expected 
to show a simple behavior with time with well-known regimes 
(e.g. Elhmaidi et al.l993; Zouari & Babiano 1994; Schumacher 
& Eckhardt 2002). The dispersion has an initial ballistic regime 
~ ^ as long as the particles lie within the viscous sub-range. For 
times larger than the Lagrangian integral scale (i.e. the time as- 
sociated with the auto-correlation function of the Lagrangian 
velocity) the relative distance between tracers increases diff'u- 
sively as in an uncorrected Brownian motion, i.e. ~ t^^^. For 



intermediate time scales the dispersion follows the Richardson 
law, ~ t^l'^ (Richardson 1926). 

Figure [121 shows the evolution between z = 30 and z = 0.5 
of the pair dispersion statistics derived with A/^ = 10^ couples of 
tracers in the same cluster run. Tracers are initially placed with 
a random sampling of the grid at z = 30 and for three diff'erent 
initial separations between the pairs: |Zo = X{t = 0)| = 36, 72 
and 288 kpc. After ^ l.SGyr the trend of the pair dispersion 
loses any dependence on the initial separation, and the time 
evolution of the pair dispersion becomes consistent with a oc 
t^^^ scaling. 

This a general behavior found in our clusters, as shown in 
the left panel of Fig [131 As soon as the cluster formation starts 
(i.e. z < 2), the pair dispersion increases following a behavior 
consistent with a ~ r^^^ scaling, although significant scatter and 
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Fig. 14. Evolution of projected positions of tracer families for run H3 (top panels), HI (central panels) and H5 (lower panels). 
Every color refers to a different shell of origin. Overlaid are the contours of the projected gas density (the first contour is for 
Pgas = ^^Pcr/fb, whcrc Pcr IS the Critical density of the Universe and ft is the baryon fraction; contours are spaced with ^/2 
intervals), the side of the image and the LOS are \5Mpc (top and central panels) and \2Mpc (lower panels), respectively. 



episodic variations across this scaling are found. In the same 
panel, we additionally show the evolution of the pair dispersion 
statistic for cluster H7, by using the implicit second order in- 
terpolation scheme introduced in Sect l3.1l (Eql4b: the reported 
trend is independent of the time interpolation scheme, and the 
scatter between the two methods is much smaller than the pair 
dispersion at all times. 

Even if the mass accretion histories of clusters and their 
velocity field may be diff'erent, the emergence of a common 
behavior suggests that the transport of tracers due to gas mo- 
tions in the ICM is similar for most of the lifetime of simulated 
clusters. In the right panel of Fig [13] we show the 3-D power 
spectrum, E{k), of the velocity field at z = for all clusters of 
the data sample. The curves of E(k) are calculated with a stan- 
dard FFT-based approach, adopting a zero-padding technique 
to account for the non-periodicity of the considered volume 
(see also Vazza et al., 2009). We also verified tat the use of 
some apodizing function, which that would help to avoid spu- 
rious eff'ects due to the sharp cut at the boundaries of the AMR 
region, provides only negligible changes in the estimate of E(k) 



(and only at scales close to the Nyquist frequency of the spec- 
tra), since the average velocity fields at the boundaries of the 
AMR region are much smaller compared to the velocity field 
inside clusters. In order to emphasize the similarity between 
our clusters, the spatial frequency k has been reseated to that of 
the the virial radius of each cluster, k{) ^ 2n/Ryir. E(k) shows 
a maximum scale in the range corresponding to ~ 1 - 2Ryir 
and a well-defined power-law behavior for nearly two orders 
of magnitude in spatial scale; the flattening found from scales 
corresponding to ~ 1 /(8A) is due to the numerical dissipation 
in the PPM scheme (e.g. Porter & Woodward 1994; Kitsionas 
et al.2009). 

The power spectrum in the velocity field is due to the com- 
bination of the turbulent motions inside the cluster volume 
with the pattern of matter accretions that surrounds the clus- 
ter volume, and with the laminar infall motions that penetrate 
inwards. The maximum coherence length of this complex ve- 
locity field is constrained by the size of the clusters themselves 
and this easily explains why the power spectrum peaks at 1-2 
Ryir. Tracer particles are frozen inside this "correlated" veloc- 
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Fig. 12. Time evolution of the pair tracers dispersion for N = 
10^ tracers in run H7. Couples with three initial separations are 
considered: |Zo = Z(0)| = 36kpc (black dots), |Zo = Z(0)| = 
12kpc (blue dots) and |Zo = Z(0)| = ISSkpc (green dots). The 
additional dashed lines show two possible oc trends to guide 
the eye. 

ity field, and the distance between two tracers is always smaller 
than the maximum scale of E(k): this explains the Richardson- 
like behavior of the pair dispersion reported in Figs fT2lfT3l The 
complex velocity field in the volume of simulated clusters al- 
lows for a fairly fast transport of particles, because large scale 
motions are strong, with particles separating by 100 - 200 
kpc distances in about IGyr. We are therefore lead to the con- 
clusion that in the analyzed set of simulated galaxy clusters the 
efficiency of particle transport is much larger than that due to 
thermal diffusion (e.g.Shtykovskiy & Gilfanov 2009), and that 
the typical scale of transport is larger than the scale of turbu- 
lent diffusion induced by central AGN activity (e.g. Rebusco et 
al.2005). Therefore the effect of large scale turbulent and mix- 
ing motions in the cluster volume is likely a key ingredient in 
any practical modeling of gas mixing processes in the ICM. 

4.2. Mixing 

The "volume-averaged" statistics examined in the previous sec- 
tions show that the transport of gas particles in the ICM is sim- 
ilar from cluster to cluster. We might however question if the 
final spatial distributions of tracers show any dependence on 
the dynamical state of the host cluster. We explored this issue 
by applying a simple initial tracers setup in a subset of repre- 
sentative objects of our sample and focused on the evolution 
of radial mixing (e.g. mixing of tracers originated at different 
distance from the center of clusters) as a function of their dy- 
namical evolution. 

A number of 10^ tracers was randomly placed within the 
cluster with a number density that follows the gas density pro- 
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Fig. 15. Map of projected gas pressure (left panels) and of pro- 
jected tracers mixing (right) panels, for run H3 (top), HI (cen- 
ter) and H5 (bottom) at z = 0. The side of the images and the 
LOS are 5Mpc; the mixing maps are produced by considering 
only cells containing more than one tracer. 



file of the clusters at z = 1 . We divided every tracer generation 
into ten families, by dividing the initial distribution to ten con- 
centric shells with an equal number of tracers and we followed 
the evolution of the system by keeping the information of the 
shell of origin for each tracer. 

In Fig [HI we show maps of the tracer evolution in clus- 
ters H3, HI and H5. These objects are good representatives 
of typical evolutive histories in our simulations: H3 does not 
experience any significant accretion of gas/DM satellites af- 
ter z < 0.2, and is fairly relaxed at z = 0; HI is interested 
by continuous accretions of gas/DM satellites and presents an 
elongated structure at the final epoch of observation; H5 expe- 
riences a major merger at z ~ 0.4, and retains a very perturbed 
dynamical state even at z = 0. 
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Fig. 16. Evolution of gas entropy profile (top panels) and of 
mean mixing (bottom panels) for HI (left) and H5 (right). 
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Fig. 17. Evolution of mean mixing profile for cluster HI 
at three redshifts, for the standard mesh refinement strategy 
(dashed lines) and for our implemented mesh refinement strat- 
egy (solid lines). 



We quantified the the degree of mixing in every cell be- 
tween tracers of "s" families and the other tracers through 



Ms = l- (\ns - Yj ^ Yj 



(5) 



where n^ is the number density of the tracers within a cell 
(at the highest resolution level) and the sum refers to all the 
species of tracers. This formula generalizes the more simple 



case of mixing between two species (e.g. Ritchie & Thomas 
2002). 

The total mixing in each cell, M, is the volume average be- 
tween all species, M - Yus^sl^s^ where A^^ is the number of 
families considered. The interpretation of this parameter is sim- 
ple: for a cell where ni ^ n2 ^ n^... ~ ^lo the diff'erent families 
are well mixed and we have M ^ 1 , while M ^ implies no 
mixing within the cell. 

Figure [15] shows the projected distributions of the mixing 
parameter M at z = for the three clusters considered above 
together with projected maps of gas pressure along the line 
of sight. We found that the diff'erences in the dynamical his- 
tory of the three objects produce diff'erent spatial distribution 
of mixing at z = 0. H3 shows a uniform pattern of mixing, with 
M ~ 0.2 - 0.3 across the whole cluster core. On the other hand 
H5 presents the most patchy mixing map with a large plume 
of the size of ~ 2Mpc in the direction of the major merger and 
a maximum mixing of M ~ 0.4 between the cores of merging 
clusters. 

In Fig [16] we compare the evolution of the gas entropy pro- 
files for HI and H5 (top panels) with that of the mixing pro- 
files (bottom panels). We found that eflftcient mixing is always 
confined in central low entropy regions. The very regular evo- 
lution of mixing in HI closely reflects the regular evolution 
of the gas entropy profile. At the final epoch of observation, a 
flat mixing profile with M 0.2 was found within the entropy 
core at r < lOOkpc, and outside this core M falls quickly to 
very low values. This is not surprising: in the idealized ICM 
represented by these simulations, convective stability to radial 
displacement is provided by the classic Schwarzschild crite- 
rion, dlog(S)/dr > (Schwarzschild 1959). A very flat entropy 
core is just marginally stable to radial perturbations, whereas 
the steep rise of the entropy profile at r > lOOkpc provides a 
very stable configuration to radial perturbations. The situation 
is more complex for the major merger system H5. The merger 
event at z ~ 0.4 - 0.5 aff'ects the entropy profile and also drives 
a more extended pattern of mixing. At the final epoch a peak of 
mixing is found at r ~ 300kpc, and significant mixing is found 
even out of r = 1 - IMpc. 

This large mixing pattern in extended merger systems con- 
firms previous works that analyzed idealized 2-body encoun- 
ters with the SPH techniques (e.g. Ritchie & Thomas 2002; 
Ascasibar & Markevitch 2006), and thus our exploratory study 
extends their results to fully cosmological Eulerian simula- 
tions. In addition, as shown by Mitchell et al.(2009), major 
mergers in Eulerian AMR codes are more efficient in mixing 
the ICM compared to mergers in SPH codes, since in SPH 
mixing during the central stages of the merger is significantly 
prevented by the suppression of instabilities by the artificial 
particles viscosity; our findings seem to support the presence 
of large eddies that provide efficient mixing during the central 
stages of cluster mergers. It is important to compare results for 
the evolution of mixing obtained through our refinement strat- 
egy (Sect j2.2l) with those obtained with a standard mesh re- 
finement strategy. For this purpose we re-simulated cluster HI, 
allowing only for mesh refinement according to the gas/over- 
density criterion and applied the same initial setup of spherical 
shells of tracers to measure the evolution of mixing at all times 
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Fig. 18. Time evolution of the projected positions for metal 
tracers injected at z = 1 from a single cD galaxy in the cen- 
ter of HI. Gas density contours are generated as in Fig [141 the 
side of the images is ~ 5 M/7c. 

(FiglTTb. Even if the overall behavior is similar, the refinement 
criterion using velocity jumps and gas/DM over-density pro- 
duces a larger mixing especially at large radii, and the progres- 
sive accretion of this gas also induces large mixing in the inner- 
most cluster regions at later redshifts. This clearly shows that 
the artificial suppression of chaotic motions due to coarse nu- 
merical resolution may reduce in a sizable way the mixing in 
the simulated ICM, even within the virial volume of clusters. 

4.3. Application to the metal mixing in the ICM 

X-ray observations proved that the ICM hosts -0.6 per cent 
of heavy elements in mass, corresponding to a mean metallic- 
ity Z ~ 0.3Zo, where Z© is the solar metallicity (e.g. Werner 
et al.2008 for a up-to-date review). The distribution of metals 
in galaxy clusters is characterized by profiles that are peaked 
toward the core of cooling flow clusters and rather flat in all 
the others (Vikhlinin et al.2005; Pratt et al.2007). Two dimen- 
sional metallicity maps from nearby clusters provided evidence 
for complex and non symmetric distributions of heavy ele- 
ments (e.g. Sanders & Fabian 2006; Finoguenov et al.2006). 
Furthermore, recent observations shed some light on the sig- 
nificant evolution with decreasing redshift of the cluster metal- 
licity (Balestra et al.2007; Maughan et al.2007). 

From the theoretical point of view, several processes can 
contribute to the observed patterns and evolution of metals in 
the ICM: galactic winds, ram pressure stripping, galaxy-galaxy 
interactions, intra-cluster supernovae, intra-cluster stellar pop- 
ulations etc (e.g. Schindler & Diaferio 2008 and references 
therein). Dealing with all these multi-physics and multi-scale 
processes is a challenging task for any cosmological simula- 



tion. As a consequence, although the overall scenario of met- 
als evolution can be reproduced by simulations with an accept- 
able agreement with observations (e.g. Cora 2006; Kapferer et 
al.2007; Tornatore et al.2007; Wiersma et al.2009), the details 
of the process are still largely uncertain. 

Although our simulations do not account for any treat- 
ment of chemical elements production and evolution, our tracer 
method can provide new hints about the ways injected metals 
are spread and mixed in the ICM during cluster evolution. 

We present here a simple set of exploratory studies, where 
initial injection profiles of metal tracers are assumed to repro- 
duce the metal release from galactic winds, and where the ad- 
vection of metal tracers is followed as a function of the evolu- 
tion of the underlying gas distribution. The DM mass resolution 
in our runs is high enough to detect single galaxy-like struc- 
tures in the DM distribution by means of the same halo-finder 
algorithm mentioned in Sect l4.2l Therefore below we will use 
the term galaxy for massive DM clumps with A/^ ~ 10^ - 10^ 
DM particles. We explore three (idealized) scenarios describ- 
ing the initial injection of "metal tracers": a) a single injection 
at z = 1 from the central cluster galaxy; b) ten injections which 
are spaced in time (in the range 1 < z < 0) from the cen- 
tral cluster galaxy; c) ten injections as in case 'b' (in the range 
1 < z < 0) from the ten most massive galaxies within the AMR 
volume. 

Following Rebusco et al.(2006), we modeled the profile of 
injected metals from a galaxy with a /^-profile 

a(r) = ao{l + (r/rjV, (6) 

where ao is the solar abundance, is the scale radius for 
each galaxy, and = 0.5. With a very simple prescription we 
adopted a constant value = 36kpc (one cell) for all galaxies 
and a normalization ao that depended on the virial DM mass 
of each cluster (with ao = I for the main cluster galaxy and 
(^o,gai Mdm,gai/Mdm,mam for the Other galaxics). 

Our choice of parameters was just a trial, designed to fall 
within the range of the parameters for galaxies studied in 
Rebusco et al.(2006). With this simple setup, we explored the 
eff'ects played by spatial transport on the mixing of metals de- 
posited in the evolving ICM. 

In Fig [18] we show the evolution of the projected positions 
of N ^ 4 • 10^ metal tracers in the cluster HI for model 'a'. 

Even if the dynamical history of this cluster at z = is 
fairly regular, the amount of mergers and accretions of matter 
experienced at z < 1 is enough to spread the metals out to 
~ SOOkpc^ 

The evolution of the number density profile of metal tracers 
from a single injection at z = 1 (model 'a') is shown in the left 
panel in Fig [19] for the same cluster. This injection scenario 
produces a very flat distribution of metals at z ~ 0.1. Even if 
the qualitative behavior is similar to that reported in Rebusco et 

^ We also notice that significant offsets (~ lOOkpc) between the po- 
sition of the DM and gas peaks are common; since X-ray observations 
are usually referred to the luminosity peak to derive profiles, we stress 
that this may have an effect because the bulk of metals is released by 
galaxies. 
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al.(2005), we notice that in our case the distance where metal 
pollution is efficient in the ICM is larger, which is consistent 
with the larger transport pattern presented in Sect l4.2[ 

The result of ten injection epochs (one every ^ 7 00 My r) 
from the central galaxy (model 'b') is shown in the central 
panel in Fig[T9j This obviously increases by ~ 10 the final 
metal content in the innermost cluster region and is found to 
produce a monotonically decreasing tracer distribution at all 
radii. 

A case of more astrophysical relevance is perhaps our 
model 'c' , which assumes multiple injections from the ten most 
massive galaxies in the AMR region (see right panel in FigfT9l). 
The corresponding spatial evolution of metal tracers in clus- 
ter HI with overlaid gas density contours, is shown in Fig|2Ql 
In this case the metal distribution is quite regular at all red- 
shifts compared to the case of injection from the central cluster 
galaxy only. 

The role played by mergers in the disruption of cool cores 
is still an object of debate (e.g. Santos et al.2009 and refer- 
ences therein), although observations support the destruction of 
cool-cores by cluster mergers (Allen et al.2001; Sanderson et 
al.2006; Rossetti & Molendi 2009). Indeed, observations sug- 
gest a dichotomy between the metallicity profiles in the cool- 
core and non-cool core clusters (e.g. De Grandi et al.2004). 
Our simulations neglect any treatment of cooling and of energy 
feedback mechanisms from central AGN, and consequently we 
cannot address this in a self-consistent way. We can speculate 
though that the most relaxed clusters in our sample are simi- 
lar to cool core clusters, while clusters with mergers are more 
similar to non-cool core clusters. 

We calculated the metallicity profiles for four clusters in 
our sample (two with a fairly relaxed dynamical evolution and 
two with a major merger in the range < z < 1) to investigate 
the eff'ect of cluster dynamics on the shape of metallicity pro- 
files. The total iron mass in the ICM is normalized so that the 
metallicity at the position of the central galaxy in the cluster 
was Z = 1 ; the gas mass was directly taken from the cells in 
the simulation. In Fig|2T]we show the profiles of the four clus- 
ters at z ^ 0.1. The profiles are computed for spherical shells 
of a width of 36kpc and they were centered on the peak of 
the X-ray bolometric cluster emission; the contribute of freshly 
(un-evolved) metal tracers from the central cluster galaxy at the 
redshift of observation was removed in post-processing. 

We found a remarkable trend: "relaxed" clusters showed 
a very peaked metallicity profile, while "merger" clusters 
showed a more complex behavior with flat profiles up to the 
outer regions. Since the injection of metal tracers in all runs 
was done in the same way and with the same duty cycles and 
no physical mechanism other than hierarchical clustering was 
at work here, this test suggests that very difl'erent metallicity 
profiles in the ICM may result from the eff'ect of diff'erent mat- 
ter accretion histories. 

Future works accounting for a realistic chemical enrich- 
ment model and cooling/feedback processes may produce 
quantitative and self-consistent predictions, yet our results 
qualitatively suggest that the observed dichotomy of metallicity 
profiles might be simply explained in terms of diff'erent cluster 
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Fig. 20. Time evolution of the projected positions for metal 
tracers injected every ^ lOOMyr from different galaxies in the 
AMR region of cluster HI . The diff'erent colors refer to the dif- 
ferent epochs of injection; gas density contours are generated 
as in Fig [141 The side of the images is ~ ISMpc. 
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(crosses) at z = 0.1. 



dynamics; tracers off'er a valuable method to investigate this 
important issue. 
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Fig. 22. Simulated iron line profile for clusters HI, H3, H4 and 
H5. The red lines show the line emissivity simulated assuming 
the tracer metallicity, the black lines are for a constant metal- 
licty and the blue lines are for an average metallicity which 
scales with radius. 



4.4. Iron line profile. 

Turbulent gas motions can lead to a Doppler shift broaden- 
ing of X-ray emitting metal lines, in excess to their intrinsic 
width. Measuring the broadening of metal lines in the future 
will likely provide the most direct evidence of subsonic turbu- 
lent motions in the ICM In a number of works (Sunyaev et 



^ Very recently Sander et al. (2009) for the first time placed direct 
limits on the turbulent broadening on the emission lines in a cool core 
galaxy cluster with XMM-Newton. The ratio between the turbulent to 
thermal energy density in the cluster core is found to be less than 13 



al.2003; Inogamov & Sunyaev 2003; Dolag et al.2005) it was 
proved that if turbulent motions are at the level of what is pro- 
duced in cosmological numerical simulations, the shape of the 
Fe XXIII line would be enough modified to allow studies of tur- 
bulence in the ICM. Here we extend these works by taking into 
account the 3-D metallicity distribution that comes out in our 
simulations adopting the 'c' model as fiducial injection model 
for metal tracers. 

We thus estimated the emission along the line of sight of 
the 6.702 keV iron line with columns through the center of the 
four galaxy clusters in our sample. For every tracer we assumed 
a line emissivity proportional to the second power of its metal 
content (which was derived as in Sect j4.3l) , and for simplic- 
ity we neglected any dependence on the gas temperature. The 
eff'ect of thermal broadening and the blending with other iron 
lines is also neglected. As a control check, we also computed 
two additional iron emissivities directly from the cells of the 
simulations: in the first case we simply assumed a fixed metal- 
licity for all cells, and in the second case we assumed a uniform 
metallicity radial profile in the form of Z/Z© oc (r/7?v/r)~^^^ (e.g. 
Vikhlinin et al.2005). In all cases, the velocity of the center of 
mass of every cluster has been removed. 

Figurj22l shows the simulated iron line profile, for clusters 
HI, H3, H4 and H5 at z = 0. For every cluster the signal was 
extracted along columns with a section of 300kpc/h and an 
energy resolution of AE = 2eV was assumed. All emissions 
were normalized to have the same total luminosity =1 within 
the column. 

For clusters without a major merger, the centroid of the iron 
line and the general shape of the line from the tracers agreed 
reasonably with those from the cells. In the case of the major 
merger system (H5) the centroids and the shapes of the lines 
were quite diff'erent. Consistently with what was found in the 
previous sections, this is because mergers create complex and 
patchy distributions of metal tracers, which cannot be fully de- 



percent within 30kpc from the cluster center, which is consistent with 
the typical turbulent to thermal energy ratios found in the centres of 
our simulated clusters, e.g. Vazza et al.(2009). 
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Fig. 23. Simulated iron line profile for cluster H5 at z = 0.5. 
The meaning of lines is the same of Figl22l 

scribed by a simple radial profile. This is further highlighted in 
FigESl which is obtained at the epoch of the major merger in 
H5. 

We conclude that the realistic modeling of the iron line 
broadening requires an appropriate modeling of the distribu- 
tion of metals in the cluster volume. 

Tracers also allow for the decomposition of the convoluted 
iron line profile in the contributions from the diff'erent genera- 
tions of metal tracers. Figure [24l shows the comparison of the 
emissivity profile for cluster HI, H3 and H5 (in this case the 
section of the column is set to 600kpc/h to also highlight the 
contribution from tracers distributed at larger distances from 
the cluster center) and the relative contributions from diff'erent 
generations of tracers. The case of the the major merger system 
highlights that the signal from the total iron line is the convo- 
lution of very diff'erent emissions coming from diff'erent gener- 
ations of metal tracers and that the earlier generations are re- 
sponsible for much broader and asymmetric components com- 
pared to the latest generations that were injected after the major 
merger. 

5. Conclusions 

In this work we studied the mixing properties of the ICM in 
simulated galaxy clusters, using tracer particles that are ad- 
vected by the gas flow during the evolution of cosmic struc- 
tures. This approach allowed us to obtain valuable information 
on the Lagrangian properties of the ICM, and allowed progress 
in our understanding of the transport processes in Eulerian sim- 
ulations of galaxy clusters. Seven galaxy clusters with total 
masses in the range M ~ 2 - 3 • lO^^M©//? and difl'erent dy- 
namical states were simulated with an updated version of the 
ENZO code (Norman et al.2007), where an additional refine- 
ment criterion was designed to reach good spatial accuracy at 



shock waves up to ~ 27?v/r from the center of galaxy clusters 
(Vazza et al.2009). In a post-processing, we injected and fol- 
lowed passive tracers, and tracked their spatial evolution with 
high time resolution. Numerical tests were presented and dis- 
cussed to find the optimal strategies to inject tracers, update 
their positions in time and sample the galaxy cluster volume 
(SectO. We then applied this tracer technique to some inter- 
esting problems concerning the mixing of the ICM. 

First, we showed that the pair dispersion statistics has a well 
constrained behavior in time: for most of the cosmic time the 
separation between close couples of tracers shows a very reg- 
ular evolution in all clusters, consistently with the basic oc 
scaling expected for a simple turbulent model (Sect j4.1l) . This 
finding is well consistent with the independent measurement 
of 3-D velocity power spectrum for the same clusters, and con- 
firms that a sizable amount of complex subsonic motions corre- 
lated on cluster scales account for 10-30 per cent of the thermal 
energy inside Ryir of our simulated clusters. 

Second, we quantified the degree of radial mixing of the 
ICM by following the evolution of tracers injected in regular 
shells within the cluster volume at z ~ 1. The radial mixing 
morphology was found to be strongly dependent on the dy- 
namical history of every cluster. Large plumes (with the size 
of ~ Mpc) where gas matter is efl&ciently mixed were found 
for a major merger cluster. On the other hand more regular and 
concentrated mixing profiles were measured for more relaxed 
clusters. 

Third, we assumed a simple model of metal injection from 
galaxies identified in the simulations around the clusters and 
studied the metal enrichment of the ICM by tracking the ad- 
vection of "metal" tracers (tracers with a diff'erent iron content 
as a function of their galaxy of origin) deposited in the cluster 
volume. In this case also, the dynamical state of galaxy clusters 
and its matter accretion history were found to aff'ect the final 
distribution of metal tracers at evolved epochs. In particular, 
the simulated clusters showed a hint of dichotomy in metal- 
licity distributions between clusters experiencing only minor 
mergers and clusters with major mergers, with the latter show- 
ing a broader distribution of metals. 

Finally, we simulated the emissivity of the Fe XXIII line 
from the center of our clusters, using the 3-D distribution of 
metals produced through our metal tracer injection (model 'c'). 
We reported significant departures from the simplified assump- 
tion of a constant or a radially symmetric metallicity profile, 
particularly for merger systems. This stresses the importance 
of studying in detail the connection between the metal enrich- 
ment history and the matter accretion history of the ICM. 



In conclusions, we believe that the coupling of high res- 
olution Eulerian numerical simulations and a tracer technique 
represents a very powerful tool to study crucial open problems 
concerning the mixing properties of the diff'use gas in galaxy 
clusters, and the injection and the advection of cosmic rays par- 
ticles in galaxy clusters, which will be subject of future works. 
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